Model vs. satellite sea ice thicknessΒΆ
In this notebook, we compare remote sensing-derived sea ice thickness from ICESat-2 and model-derived sea ice thickness from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS, Zhang and Rothrock, 2003)
import xarray as xr # For working with gridded climate data
import pandas as pd
import numpy as np
import itertools
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.misc_utils import restrictRegionally # Region restriction
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting
# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook
# Remove warnings to improve display
import warnings
#warnings.filterwarnings('ignore')
1) Read in dataΒΆ
Below, weβll read in the data and restrict it to the Inner Arctic region. For the PIOMAS sea ice thickness, weβll mask out the pole hole and set 0 thickness to nan to match the ICESat2 data; the PIOMAS data sets anywhere with no sea ice to nan, whereas PIOMAS sets anywhere with no sea ice to 0.
book_ds = read_book_data() # Read/download the data
book_ds = book_ds[["ice_thickness_smoothed","piomas_ice_thickness"]] # Grab data variables of interest
book_ds = restrictRegionally(book_ds, regionKeyList=[10,11,12,13,15]) # Restrict data to the Inner Arctic
book_ds["piomas_ice_thickness"] = book_ds["piomas_ice_thickness"].where(book_ds["piomas_ice_thickness"]!=0, np.nan) # Set 0 --> nan
book_ds = book_ds.where(book_ds.latitude<=88, np.nan) # Mask out pole hole
years = [2018,2019,2020] # Years over which to perform analysis
Regions selected: Inner Arctic
2) Map mean winter differencesΒΆ
Here, weβll compute winter means for PIOMAS and ICESat-2 sea ice thickness. Then, weβll use the interactiveArcticMaps function to generate maps of the winter means, along with a comparison map showing the gridcell differences.
is2_winter_means = compute_gridcell_winter_means(book_ds["ice_thickness_smoothed"], years=years, start_month="Nov")
pio_winter_means = compute_gridcell_winter_means(book_ds["piomas_ice_thickness"], years=years, start_month="Nov")
for i in range(len(years)):
data1 = is2_winter_means.isel(time=i)
data2 = pio_winter_means.isel(time=i)
pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=data1.time.values.item())
pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title=data2.time.values.item())
diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference")
data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
display(data_pl)
print("\n")